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Abstract 

Many body trial wave functions are the key ingredient for accurate Quantum Monte Carlo estimates of total electronic 
energies in many electron systems. In the Coupled Electron- Ion Monte Carlo method, the accuracy of the trial function 
must be conjugated with the efficiency of its evaluation. We report recent progress in trial wave functions for metallic 
hydrogen implemented in the Coupled Electron-Ion Monte Carlo method. We describe and characterize several types 
of trial functions of increasing complexity in the range of the coupling parameter 1.0 < r s < 1.55. We report wave 
function comparisons for disordered protonic configurations and preliminary results for thermal averages. 
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1. Introduction 

Modern ab initio simulation methods for sys- 
tems of electrons and nuclei mostly rely on Density 
Function Theory (DFT) for computing the elec- 
tronic forces acting on the nuclei, and on Molecular 
Dynamics (MD) techniques to follow the real-time 
evolution of the nuclei. Despite recent progress, 
DFT suffers from well-known limitations [1,2]. As a 
consequence, current ab initio predictions of met- 
allization transitions at high pressures, or even the 
prediction of structural phase transitions, are often 
only qualitative. Hydrogen is an extreme case [31415] . 
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but even in silicon, the diamond//?-tin transition 
pressure and the melting temperature are seriously 
underestimated [6] . 

An alternative route to the ground-state prop- 
erties of a many-electrons system is the Quantum 
Monte Carlo method (QMC) 02]. In QMC, a many- 
body trial wave function for the electrons is as- 
sumed and the electronic properties are computed 
by Monte Carlo methods. Bosonic details of the 
trial wave functions are automatically optimized by 
projecting the trial state onto the ground state with 
the same nodal structure. Hence for fermions, QMC 
is a variational method with respect to the nodes 
of the trial wave function and a systematic, often 
unknown, error remains |7l2j . Over the years, the 
level of accuracy of the fixed-node approximation 
has been improved [819|10|llj such that, in most 
cases, fixed-node QMC methods have proven to be 
more accurate than DFT-based methods, on one 
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side, and less computationally demanding than cor- 
related quantum-chemistry strategies (such as cou- 
pled cluster method) [2] on the other side. Recently 
we have developed an ab-initio simulation method, 
the Coupled Electron-Ion Monte Carlo (CEIMC) 
method, based entirely on Monte Carlo algorithms, 
both for solving the electronic problem and for 
sampling the ionic configuration space in the Born- 
Oppenheimer approximation [12]. A Metropolis 
Monte Carlo simulation of the ionic degrees of free- 
dom (represented either by classical point particles 
or by path integrals) at fixed temperature is per- 
formed based on the electronic energies computed 
during independent ground state Quantum Monte 
Carlo calculations. Application of CEIMC has been 
limited so far to high pressure hydrogen for several 
reasons: a) hydrogen is the simplest element of the 
periodic table, and the easiest to cope with since the 
absence of the additional separation of energy scales 
between core and valence electrons as in heavier 
elements; b) it is an important element since most 
of the matter in the universe consists of hydrogen; 
c) its phase diagram at high pressure in the inter- 
esting region where the metallization occurs is still 
largely unknown because present experiments are 
not able to reach the relevant pressures. We have 
investigated the very high pressure regime where all 
molecules are dissociated and the system is a plasma 
of fully ionized protons and electrons [13] , and we 
have studied the pressure-induced molecular disso- 
ciation transition in the liquid phase [T4], In both 
studies the CEIMC results were not in agreement 
with previous Car-Parrinello Molecular Dynamics 
(CPMD) calculations [Hill]. While we have evi- 
dence now that the discrepancy in the fully ionized 
case is removed by taking a more accurate trial wave 
function in the QMC, in the second study more 
accurate CEIMC calculations predict a continuous 
molecular dissociation with increasing pressure at 
variance with CPMD where a first order molecular 
dissociation transition was observed by increasing 
pressure at constant temperature[16] . More re- 
cently, using constant volume Born-Oppcnhcimcr 
Molecular Dynamics rather then constant pressure 
CPMD, a continuous dissociation transition has 
been reported from DFT-GGA studies [P7] . 

In the present paper we discuss in some details the 
various trial wave functions we have implemented 
for hydrogen and we compare their accuracy at var- 
ious densities. We will not review the details of the 
CEIMC method which have been described at length 
in reference [12]. We just mention that in metals 



huge finite size effects, caused by the discrete na- 
ture of the reciprocal space, can be alleviated by 
averaging electronic properties over the overall un- 
determined phase of the many-body wave function 
(Twist Average Boundary Conditions). Results re- 
ported in this work are obtained by this method 
(see also ref.[TH] for recent development). Section [2] 
will be devoted to describing the different trial wave 
functions and some details on their efficient imple- 
mentation. In Section [3] we will report numerical 
comparisons among the various trial functions. Fi- 
nally, in Section 0] we collect our conclusions and 
perspectives. 

2. Trial Wave Functions for Hydrogen 

The trial wave functions we have adopted for hy- 
drogen are of the simple Slater-Jastrow form. We 
have considered spin unpolarized hydrogen only. For 
each spin state, a single determinant of one electron 
orbitals 0fe(r) is used to account for the fermionic 
symmetry of the many-body wave function. A Jas- 
trow factor e~ u is then added to account for (at 
least) two body correlations directly into the trial 
wave function 
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and U = 2j<i u iji r ij)- An accurate and parame- 
ter free Jastrow factor can be obtained within the 



RPA uf i PA {i 



e iV) [19j - This simple form 



satisfies the correct cusp conditions at short particle 
separations and the right plasmon behavior (screen- 
ing) at large distances. It was shown [213] to pro- 
vide good energies for hydrogen even at intermedi- 
ate densities if supplemented by gaussian functions 
Uij{r) = u R j PA (r) — ctij e~ r , with the variational 
parameters onj , Wij . The additional terms preserve 
the short- and long-distance behavior of the RPA 
function and correct for possible inaccuracies at in- 
termediate distances. However, they introduce four 
variational parameters to be numerically optimized, 



namely a 



We have adopted the RPA 
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form of the two body correlation in all different trial 
function discussed below. 

In Ref. [10] we have described a procedure to it- 
eratively improve any initial trial function based on 
the Feynman-Kac formula. If a determinant of single 
electron orbitals is assumed as an initial ansatz, the 
first iteration generates a bosonic (symmetric) two- 
body correlation function (Jastrow) while the next 
iteration naturally provides the backflow transfor- 
mation of the orbitals and a three-body bosonic cor- 
relation term. Unfortunately, this is a formal theory 
which cannot provide, in general, analytical expres- 
sions for the various terms. Nonetheless the general 
structure is illuminating in searching for improve- 
ments. 

In the first implementation of CEIMC [2TT22] m 
single electron orbitals consisting of a linear com- 
bination of few optimizable guassian orbitals cen- 
tered on each molecule was used in order to simulate 
the insulating phases of molecular hydrogen. Opti- 
mization of the variational parameters, in number 
proportional to the number of electrons in the sys- 
tem, needed to be performed at each ionic configura- 
tion and was a major bottleneck for the efficiency of 
the method. Subsequently, we have developed trial 
functions with a very limited number of variational 
parameters (even zero when possible) and therefore 
reduced the complexity of the optimization step in 
the electronic calculation (or to reduce to a linear 
optimization in the case of DFT orbitals). 

2.1. The Metallic Wave Function 

At a very high pressure, beyond metallization and 
molecular dissociation, the electron liquid is a good 
Fermi liquid and correlation effects, with protons 
and among electrons, can be treated as perturba- 
tions. In this case it is natural to assume a determi- 
nant of free electron states (plane waves) as an initial 
ansatz. As stated above, the first iteration provides 
the two-body Jastrow factor for which we adopted 
the modified RPA form discussed above. Although 
this form of the trial function has been extensively 
used in the early QMC study of hydrogen [20], it 
is worthwhile to emphasize that the proton coor- 
dinates do not appear in the fermionic part of the 
trial function, that is in its nodal structure, which is 
ultimately the limit of the accuracy of QMC. This 
limitation is overcome at the next iteration in the 
Feynman-Kac procedure, which suggests the back- 
flow transformation of the orbitals and a three-body 



correlation factor. The form of the backflow trans- 
formation is 

x 4 = n + (M) r U' ( 3 ) 

j 

where r} a p are the electron-electron and electron- 
proton backflow functions that must be parameter- 
ized. When the single body-orbitals in the determi- 
nants are expressed in terms of the quasi-particle 
coordinates Xj, the nodal surfaces of the trial wave 
function become explicitly dependent on the proton 
positions, a crucial characteristic for inhomogeneous 
electron systems which will provide a more accurate 
energy even at the RQMC level. Similar to the case 
of the homogeneous electron gas [5] , the backflow and 
three-body functions were at first parametrized as 
gaussians[22]. This trial function has a total of 10 
free parameters to be variationally optimized and 
has been used in a first CEIMC study of the melting 
transition of the proton crystal in hydrogen at r s = 
1 [22] . Subsequently, we were able to derive approxi- 
mate analytical expressions for the backflow and the 
three-body functions, as well as for the two-body 
correlation factor, in the Bohm-Pines collective co- 
ordinates approach [lOj . This form is particularly 
suitable for the CEIMC because it is parameter-free. 
At the same time, it provides comparable accuracy 
to the numerically optimized wave function, both in 
the crystal configuration and for disordered protons. 
Explicit forms of the various terms can be found 
in the appendix of Ref. [10], With this kind of wave 
function we have investigated the melting at three 
densities (r s = 0.8,1.0,1.2) including quantum ef- 
fects for protons [13| 12] , 

2.2. Band- structure-based Wave Functions 
(IPP/LDA ) 

The metallic wave function is expected to pro- 
vide an accurate description of the electronic ground 
state at high density, well beyond molecular disso- 
ciation and metallization. On the other hand we ex- 
pect it to be a poor representation of the true ground 
state at lower densities near the molecular dissocia- 
tion. At lower densities, molecules appear and plane- 
wave single-electron orbitals (although in terms of 
backflow coordinates) are certainly not a good rep- 
resentation. Natoli et aL [23124] have previously used 
a Slater determinant of Kohn-Sham self consistent 
orbitals to study the solid phases of atomic hydrogen 
at r s = 1.31 and T = 0, and of molecular hydrogen 
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at lower densities. They have found a typical energy 
gain of 0.5eV/electron by replacing the plane- wave 
with the self consistent orbitals in the Slater deter- 
minant. Here we have implemented similar ideas. 

The single-particle orbitals, {4> n }, that comprise 
the Slater determinant are computed on-the-fly dur- 
ing the CEIMC calculation as the eigenstates of 
some single-particle Hamiltonian 



-G 2 < E c 



(10) 



In the plane-wave basis, the orbitals can be ex- 
pressed as 



„ k (r) = ^£C nkG e*( k + G K 

V v G 



h(f> n (r) = e n (j) n (r) , 



(4) 



where the N/2 orbitals with the lowest eigenvalue, 
e n , are selected to fill the determinants. The single- 
particle Hamiltonian describes electron-nuclear in- 
teractions and approximates electron-electron inter- 
actions through an effective potential. In Hartree 
atomic units 



(11) 



The problem of filling the Slater determinant re- 
duces to a problem of finding the plane-wave co- 
efficients, C„kG, of the Bloch states and evaluat- 
ing the resulting determinant for the QMC elec- 
tronic configuration. The plane-wave orbital coeffi- 
cients are found by solving the eigenvalue problem 
of the single-particle Hamiltonian 



h= -iv 2 + ^ ff (r; S) 



(5) 



G' 

where 



5G 



Be ) C'rikG' = ZnkCnkG, 
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The effective potential, V c g(r;S), depends on the 
nuclear coordinates, S, as parameters and is, by con- 
struction, invariant under translation of the electron 
coordinate from one simulation cell to another: 



B, 



G 



B G - > = -||k + G| 2 5 G , G , 



(13) 



V eS (r + L; S) = V cS (r; S) . 



(6) 



i J V cS (r;S)e< G - G >d 3 r 



The many-body Slater- Jastrow wave function 
should obey twisted boundary conditions [25 . The 
effect is that translating any single-electron coordi- 
nate by a lattice vector twists the overall phase of 
the wave function 



* (ri 



r % + L, . . . , r N ) = e 1 "* (n, . . . , n, . . . , rjv) (7) 



which is achieved by enforcing Bloch's theorem on 
the single-particle orbitals for the first-Brillioun 
zone wavevector corresponding to the twist angle 



0nk (r) = Unk (r) e lk r ;ke 1BZ. 



(8) 



The Bloch functions, U n k, are then cell-periodic 
functions, which are particularly well-represented 
in a plane- wave basis set, with basis functions of 
the form 



iG.: 



(9) 



The basis functions are orthogonal in G and nor- 
malized to 1 through the simulation cell of volume 
V. The G vectors are vectors of the reciprocal-space 
lattice, defined such that G.L = 2rmr;m 6 Z. The 
quality of the plane- wave basis is controlled by a sin- 
gle parameter, an energy cutoff, which determines 
that largest wavevector that is present in the set: 



From this Fourier transform, it is clear that if the 
simulation cell is centrosymmetric then the poten- 
tial is real and inversion symmetric in G — G'. How- 
ever, this and other symmetry considerations, such 
as the existence of point-group operations that leave 
the single-particle Hamiltonian invariant, are gen- 
erally not relevant in a CEIMC simulation in which 
the nuclear coordinates are typically not ordered. 
One relevant symmetry, time-reversal symmetry for 
k — > — k, implies that only half of the Brillioun zone 
must be explicitly sampled. Further, the eigenvec- 
tors corresponding to k = are real. 

Note that the effective potential in reciprocal 
space depends on G — G' only, and can therefore 
be stored as a one-dimensional vector of size 8Nq 
in place of an array of size Nq. This property is due 
to the locality of the potential in real space. One 
method of solving the eigenvalue problem would be 
to use standard numerical matrix diagonalization 
techniques. However, such approaches are much too 
slow when the number of basis functions signifi- 
cantly exceeds the number of required eigenstates, 
as in the case of a plane- wave basis. Instead, we use 
an iterative, conjugate-gradient band-by-band mini- 
mization scheme [2"S] that is particularly effecient for 
these classes of problems. The method employs the 
variational principle to minimize residuals and the 
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Gram-Schmidt scheme to preserve orthogonaliza- 
tion of the eigenstates. In our studies, we typically 
use one of two choices for the effective potential: 

- V e s = V e - n , the bare electron- nuclear interaction. 
We refer to this method as the independent par- 
ticle potential (IPP) because electron-electron in- 
teractions and nuclear screening are absent from 
the orbitals. 

- V e ff = Vlda, the Kohn-Sham effective potential 
within the local density approximation (LDA) . 

Further complexities may be introduced to all of 
these wave functions through backflow transforma- 
tions and the Jastrow factor. 

2.3. IPP and LDA Orbitals 

For hydrogen, the IPP orbitals are generated with 
V c f{ from the bare electron-proton Coulomb poten- 
tial such that the potential term from Equation [T"4l 
becomes 

V ep (G - G') = is G _ G ,v c (|G - G'|) , (14) 
where Sg is the simulation-cell structure factor 
S G =5>- iG - R *, (15) 

with Ri a proton coordinate, and v c is the Fourier 
transform of the Coulomb potential 



v c (G) 



47T 
Q2- 



(16) 



The diverging G = G' components are equated 
to zero, which corresponds to shifting the arbitrary 
zero of electrostatic potential. 

DFT-LDA orbitals are generated from the Kohn- 
Sham effective potential with the local density ap- 
proximation. The effective potential becomes a sum 
of three terms 

Kff (r) = Vep (r; S) + V Ka (r) + V^ A (r) . (17) 

Both the Hartree potential and the exchange- 
correlation potential depend on the electron density 



(18) 



which in turn depends on the eigenfunctions of the 
Hamiltonian. These potentials must therefore be de- 
termined self-consistently, through an iterative pro- 
cedure, such that the Hamiltonian yields orbitals 



corresponding to a density equal to that used to gen- 
erate the potentials. The local density approxima- 
tion that we use is the Perdew-Zunger[27] parame- 
terization of Ceperley- Alder [28J electron-gas data. 

For hydrogen, both the LDA and IPP wave func- 
tions are eigenstates of a Hamiltonian which con- 
tains a bare Coulomb interaction between electrons 
and protons. The singularity in the potential results 
in a derivative cusp in <j> (r) when r = R/. The pres- 
ence of this cusp with orbitals that properly sat- 
isfy Kato's cusp condition is important for obtain- 
ing good energies or short projection times for QMC 
algorithms. The electron-proton cusp condition is: 



d<f> (r) 



Or 



»(r = Rj) 



(19) 



Representing this cusp in reciprocal space is chal- 
lenging due to the slow algebraic decay of G~ 4 . For 
this reason, we implement a cusp-removal method 
by dividing the orbitals by a function that satisfies 
the cusp condition exactly (in this case, the RPA ep 
Jastrow function discussed earlier), 



</>nk (r) 



<^nk (r) 



(20) 



before reverse Fourier transforming to retrieve the 
C„kG coefficients used to build the Slater determi- 
nant. The proper electron- proton cusp is later rein- 
troduced exactly in real space using the same RPA 
Jastrow function. This procedure greatly enhances 
the convergence of the Slater- Jastrow wave function 
with respect to the size of the plane-wave basis set 
used to represent the orbitals, as demonstrated by 
Figure Q] 

2.4. Backflow Transformations of IPP /LDA 
Orbitals 

Similarly to the metallic wave function case, we 
can think of applying the Feynman-Kac iteration 
procedure to generate a backflow transformation 
even for IPP/LDA orbitals. This transformation 
introduces correlation effects into the Slater deter- 
minant part of a Slater- Jastrow wave function, with 
the advantage that modifications of the nodal sur- 
face are possible and the energy is improved. The 
specific form for the transformation is unknown, 
but as a first ansatz we can use the same analyti- 
cal expressions we have developed in the metallic 
case. As we will show in the next section, backflow 
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Fig. 1. VMC total energy for a fixed proton configuration 
and a Slater- Jastrow trial wave function with IPP orbitals. 
The VMC energy converges slowly with respect to the or- 
bital plane-wave basis cutoff, but this is greatly improved 
by removing the inaccurate ep cusp and correcting in real 
space with an analytic Jastrow function. 

is found to improve the total energy and the vari- 
ance of the total energy, both indicating that it is a 
better representation of the ground state of the sys- 
tem. Only the electron-electron backfiow is applied, 
while the electron-proton term is found to reduce 
the quality of the trial function. 

3. Comparisons of Wave Functions 

3.1. Fixed protonic configurations 

In this section we consider fixed protonic config- 
urations and we compare the quality of the vari- 
ous wave functions at several densities. We consider 
first protons in crystal structures at r s = 1.0 and 
r s = 1.4. Results for BCC hydrogen at r s = 1.31 
obtained from various improvements of the metal- 
lic wave function are reported in Table III of Ref. 
[TO] where they are also compared to the results ob- 
tained with self-consistent Kohn-Sham orbitals [23] ■ 
There we have shown that the quality of the analyt- 
ical form of the metallic wave function is superior to 
its numerically optimized version and comparable to 
that of the LDA orbitals for hydrogen in the BCC 
structure and for various system sizes. Our present 
implementation of LDA orbitals provides results in 
agreement with previous estimates |29). 

In Tables [JJ and [5] we report QMC energies for 
hydrogen in several crystal structures and for var- 
ious system sizes at two densities corresponding to 
r s = 1 and r s — 1.4 respectively. A complete study 



Table f 

Energy and variance of hydrogen in various structures with 
different trial functions at r s = f .0. All results are obtained 
averaging over a 6x6x6 fixed grid of twist angles. E v and 
<r^ represent, respectively, energy and variance at the varia- 
tional level while E r and a^. are the energy and the mixed 
estimator for the variance obtained with RQMC. Units are 
Hartree/atom. 





WFS 


E v al E r cr 2 r 


BCC 

(N p = 54) 


Met 

IPP 

LDA 
IPP+BF 
LDA+BF 


-0.36931(1) 0.0279(2) -0.3721(1) 0.0182(7) 
-0.3681(3) 0.0765(3) 
-0.3681(2) 0.0765(2) 
-0.3705(1) 0.0359(1) 

-0.36805(2) 0.04636(4) -0.3705(1) 0.0357(1) 


FCC 
(N p = 32) 


Met 

IPP 

LDA 
IPP+BF 
LDA+BF 


-0.3792(1) 0.01543(4) 
-0.3756(2) 0.0756(2) 
-0.3757(2) 0.0753(2) 
-0.3779(1) 0.03550(9) 
-0.3779(1) 0.0351(1) 


DIAM 

(N p = 64) 


Met 

IPP 

LDA 
IPP+BF 
LDA+BF 


-0.3477(2) 0.0268(1) 
-0.3621(2) 0.0830(4) 
-0.3613(2) 0.0823(3) 
-0.3637(1) 0.0404(1) 
-0.3635(1) 0.0406(1) 


DIAM 

(N p = 8) 


Met 

IPP 

LDA 
IPP+BF 
LDA+BF 


-0.41060(4) 0.02136(4) -0.41368(6) 0.01032(2) 
-0.40198(8) 0.0863(1) -0.4094(1) 0.04342(8) 
-0.40206(8) 0.0865(1) -0.4098(1) 0.04356(6) 
-0.40632(6) 0.04958(8) -0.41070(6) 0.02382(4) 
-0.40638(6) 0.04958(8) -0.4107(1) 0.02382(6) 



of the size dependence and the relative stability of 
those structure is not our concern here and will be 
reported elsewhere [29] . We observe that LDA always 
provides a small or negligible improvement over IPP, 
while IPP is significantly cheaper through the lack 
of the self-consistent requirement. Comparing var- 
ious structures and system sizes, we observe that 
the best wave function depends on the structure: for 
BCC, FCC structures and the diamond structure 
with N=8, the metallic wave function is superior 
to the others. The opposite is true for the diamond 
structure with N=64, where IPP and LDA provide 
lower energies at all densities. We also observe that 
the ordering of the wave functions does not appear to 
depend on density, at least in the limited range inves- 
tigated. Note that r s = 1.31 corresponds to the den- 
sity predicted by ground state QMC calculations [213] 
for the molecular dissociation to occur. Another im- 
portant test is the effect of backfiow (BF) on the 
band orbitals. We see that both variation and rep- 



G 



Table 2 

Energy and variance of hydrogen in various structures with 
different trial functions at r s = 1.4. All results are obtained 
averaging over a 6x6x6 fixed grid of twist angles. E v and 
<rj represent, respectively, energy and variance at the varia- 
tional level while E r and tr^ are the energy and the mixed 
estimator for the variance obtained with RQMC. Units are 
Hartree/atom. 



-0.35 - 





WFS 


E v al E r al 


BCC 
(N p = 54) 


Met 
IPP 
LDA 
IPP+BF 


-0.5203(2) 0.036(2) -0.5224(1) 0.01008(4) 
-0.5139(1) 0.0887(2) -0.5209(2) 0.0284(1) 
-0.5145(1) 0.0869(3) -0.5210(2) 0.0286(1) 
-0.5228(1) 0.01413(7) 


FCC 

(N p = 32) 


Met 
IPP 
LDA 
IPP+BF 


-0.5272(1) 0.00872(3) 
-0.5210(1) 0.0835(3) -0.5256(1) 0.0276(1) 
-0.5212(1) 0.0828(3) -0.5259(1) 0.02724(9) 

-0.5280(1) 0.01352(5) 


DIAM 

[N p = 64) 


Met 
IPP 
LDA 
IPP+BF 


-0.5168(1) 0.01656(8) 
-0.5189(2) 0.1027(8) -0.5321(5) 0.0339(3) 
-0.5323(1) 0.0331(2) 
-0.5346(1) 0.01740(7) 



tation energies are slightly improved and the vari- 
ational variance is halved by the backflow, which 
means a net improvement of the trial function and 
the need of a shorter projection in imaginary time to 
reach the ground state. The high level of accuracy 
observed for the metallic wave function in crystal 
structures induced us to perform a detailed study of 
liquid hydrogen at finite temperature [13] ■ 

Next, we consider how the various wave functions 
perform on disordered protonic configurations rep- 
resentative of atomic hydrogen in the liquid state in 
the range of densities corresponding to the interval 
r s S [1.0, 1.55]. As before, all results reported here 
are averaged over a 6x6x6 fixed grid in the twist 
space. At r s = 1 we compare the Metallic wave func- 
tion with the LDABF wave function, while at lower 
densities we report data for the IPP, LDA and LD- 
ABF wave functions. 

In the left panel of Figure Owe display VMC and 
fully converged RQMC energies for 18 protonic con- 
figurations at r s = 1.0 obtained with the metalllic 
and the LDABF wave functions. Configuration is 
a 32 protons warm crystal (FCC) near melting, con- 
figuration 1 is the perfect BCC crystal with 54 pro- 
tons, configurations 2 to 12 are statistically indepen- 
dent configurations of 54 protons obtained during a 
CEIMC run at T=2000K performed with the LD- 
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Fig. 2. Total energy (left panel) and quality parameter (right 
panel) for a number of static proton configurations as ob- 
tained with the metallic and the LDABF trial functions at 
r s = 1.0. TABC with a 6x6x6 fixed grid in the twist space is 
performed. Energies are in h/atom. The quality parameter 
is defined in the text. 

ABF trial function, while the remaining 5 configura- 
tions have been obtained during a CEIMC run at the 
same temperature performed with the metallic trial 
function. Several interesting facts can be inferred 
from this figure. With the noticeable exception of 
the perfect BCC crystal, energies from LDABF wave 
function are always lower than energies from the 
metallic wave function. In particular, the fully con- 
verged RQMC energies from the metallic wave func- 
tions are above the VMC energies from LDABF trial 
function. This implies that changing the form of the 
nodes provides more energy than fully projecting 
the initial state. Why the excellent quality of the 
metallic wave function observed in perfect crystals 
is deteriorated by disorder remains unclear to us, 
but as a matter of fact, it appears that LDA nodes 
supplemented by e-e backflow perform much better 
both in the liquid state and in the crystal state with 
thermal fluctuations. Another interesting observa- 
tion concerns the dispersion of the energies over a 
set of configurations. The dispersion of the energy, 
measured by its standard deviation with respect to 
the mean value, is a measure of the roughness of the 
BO energy surface, the key quantity in determining 
the structure and the thermal properties of hydro- 
gen. The dispersions over the first 11 liquid config- 
urations (from conf. 2 to conf. 13 generated during 
a run with the LDABF trial function) is reported 
in Table [31 Dispersion of the RQMC energies with 
the metallic trial function is roughly twice as large 
than the dispersion from the LDABF trial function 
which indicates as the BO surface with LDABF is 
smoother than the other, providing a less structured 
liquid state and a lower melting temperature of the 
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Fig. 3. Total energy (left panel) and quality parameter (right 
panel) for a number of static proton configurations as ob- 
tained with the metallic and the LDABF trial functions at 
r s = 1.31. TABC with a 6x6x6 fixed grid in the twist space 
is performed. Energies are in h/atom. In the right panel 
open symbol represent VMC energies for IPP (circles), LDA 
(squares) and LDABF (triangles), respectively. RQMC en- 
ergies for the same trial functions are represented bv closed 
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Fig. 4. Total energy (left panel) and quality parameter (right 
panel) for a number of static proton configurations as ob- 
tained with the metallic and the LDABF trial functions at 
r s = 1.40. TABC with a 6x6x6 fixed grid in the twist space 
is performed. Energies are in h/atom. In the right panel 
open symbol represent VMC energies for IPP (circles), LDA 
(squares) and LDABF (triangles), respectively. RQMC en- 
ergies for the same trial functions are represented by closed 
symbols. 

proton crystal. Moreover, it is interesting to compare 
the dispersion of the VMC and the RQMC energies 
for a given trial function and a given set of config- 
urations. For the same set of liquid configurations 
and for the LDABF wave function, the VMC and 
RQMC dispersions differ by 0.09mH/at correspond- 
ing to roughly 30K, a tiny difference in the rough- 
ness of the BO energy surface. This suggest that at 
r s = 1.0 the VMC BO energy surface is likely to be 
accurate enough. In the right panel of Figure [3 we 
report for the same static configurations the qual- 
ity parameter for the two wave functions considered. 
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Fig. 5. Total energy (left panel) and quality parameter (right 
panel) for a number of static proton configurations as ob- 
tained with the metallic and the LDABF trial functions at 
r s = 1.55. TABC with a 6x6x6 fixed grid in the twist space 
is performed. Energies are in h/atom. In the right panel 
open symbol represent VMC energies for IPP (circles), LDA 
(squares) and LDABF (triangles), respectively. RQMC en- 
ergies for the same trial functions are represented by closed 
symbols. 

Table 3 

Dispersion of the energy over representative static config- 
urations at various densities. All results are obtained av- 
eraging over a 6x6x6 fixed grid of twist angles. Units are 
mHartree / atom. 





VMC 


r s 


Met 


IPP LDA LDABF 


1.0 


1.666(1) 


0.820(1) 


1.31 




0.93227(2) 0.66897(2) 0.53286(1) 


1.40 




2.10255(4) 1.01702(2) 0.74335(1) 


1.55 




2.59550(5) 2.9005(1) 1.62655(3) 




RQMC 


r s 


Met 


IPP LDA LDABF 


1.0 


1.380(3) 


0.730(2) 


1.31 




0.81059(3) 0.46996(2) 0.48619(1) 


1.40 




1.60822(5) 0.76393(3) 0.75089(1) 


1.55 




2.15821(6) 1.8410(1) 1.74226(5) 



The quality parameter of a trial function is defined 
as the negative logarithm of the overlap of the trial 
state onto its fully projected state. It is easy to prove 
that it reduces to the integral over positive imagi- 
nary time of the difference between the energy and 
its extrapolation at infinite time. The smaller the 
quality parameter the better the trial function is. We 
observe that, with the exception of the BCC crystal, 
the metallic wave function has larger values, which 
means that it is less accurate than the LDABF wave 
function. Note also, how the quality of the LDABF 
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wave function is uniform (at fixed number of parti- 
cles) through the perfect crystal and the disordered 
configurations, no matter how these configurations 
have been generated. This is an important require- 
ment to accurately predict phase transitions. On the 
other hand, the quality of the metallic wave func- 
tion on the 5 liquid configurations generated with 
this wave function is higher than on the remaining 
f liquid configurations generated in a LDABF run. 
A good trial function should have a uniform quality 
throughout the entire proton configurational space 
in order to provide an unbiased sampling. 

A similar analysis has been performed at r s = 
1.31, 1.40 and 1.55 considering 5 liquid configura- 
tions generated during a CE1MC run at T = 2000X 
with the LDABF trial function. Results are dis- 
played in Figures [3J 2] and [5] respectively. Since the 
metallic wave function is certainly not accurate at 
this density, we compare the 1PP, LDA and LDABF 
wave functions only. At all densities LDABF trial 
function provides lower energies (its VMC energy 
is very close the RQMC-IPP energy). Moreover its 
quality parameter is smaller than for the other wave 
functions because the backflow provides a consid- 
erable improvement already at the VMC level and 
the gap between the VMC and the RQMC energies 
is smaller than for the other cases. Also the quality 
parameter appears to be more uniform over the 
(small number of) configuration examined. Qual- 
ity parameters of IPP and LDA trial functions are 
within error bars to each other except at r s = 1.55 
where the IPP one is quite smaller, and close to 
the LDABF values, than the LDA values. This is 
caused by an larger energy gap between VMC and 
RQMC energies and a slower convergence in imagi- 
nary time for LDA trial function than for the other 
trial functions. 

Energy dispersions are collected in Table El We 
observe a large dependence of the dispersion (be- 
tween 2 and 3 times) on the trial function and, for 
IPP and LDA, on the level of optimization of its 
bosonic part. As mentioned above, this fact would 
imply a large sensitivity of thermal properties to the 
type of trial function and on the QMC method ex- 
ploited to obtain the BO energy surface. These con- 
clusions are in agreement with our finding that the 
nature of the dissociation process changes from the 
VMC to the RQMC energy surface [14]. In the LD- 
ABF case the difference in dispersion between VMC 
and RQMC energies corresponds to only few tens of 
degrees Kelvin, probably a negligible effect in most 
interesting cases. On the other hand VMC remains 




r/a r/a 



Fig. 6. Left panel: r s = 1, T = 1500^, N p = 54. Proton-pro- 
ton pair correlation functions as obtained with LDABF and 
metallic wave functions. Results are obtained with TABC 
by using twist sampling around a 4x4x4 grid. CPMD data 
from ref. |15| are also represented by a thick dashed line. 
Right panel: r s = 1.4, T = 2000.fi:, N p = 54. Proton-proton 
pair correlation functions as obtained with IPP and LDABF 
wave functions. Results are obtained with TABC by using 
a 6x6x6 fixed grid (IPP) and by twist sampling around a 
4x4x4 grid (LDABF). 

much more efficient in terms of computational re- 
sources required. 

3.2. Liquid- State Simulations 

After the validation of LDABF trial function of 
the previous section, we report here results for the 
liquid structure of hydrogen at the same densities. 
We have simulated systems of 54 protons. The 
TABC is performed here using twist sampling [TS] 
around the nodes of a 4x4x4 grid in twist space at 
each protonic step. In the left panel of Figure [5] we 
report a comparison of proton-proton pair corre- 
lation functions g pp (r) at r s = 1 and T = 1500X 
as obtained from the metallic and LDABF trial 
functions at the VMC level. As expected from the 
results of the previous section, we observe consider- 
ably more structure with the metallic trial function 
than with the LDABF one, which indeed would 
correspond to having an effective lower tempera- 
ture. On the same figure we report results from a 
CPMD simulation [T5] performed within the LDA 
approximation. The agreement between CPMD 
data and our present CEIMC data from LDABF 
trial function is striking and somehow unexpected. 
Indeed our representation of the electronic ground 
state is much more accurate than the simpler LDA 
one. Also the finite size effect in the CPMD calcu- 
lation was addressed only partially by using only 
closed shells systems at the T point. Nonetheless 
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the observed agreement testify that the structure 
of the proton liquid is not very sensitive to details 
of the ground state representation, at least at such 
high density. Finally, in the right panel of Figure [6] 
we report preliminary data for g pp (r) of 54 protons 
at r s = 1.4 and T = 2QQQK. We compare IPP and 
LDABF trial functions at the VMC and RQMC 
levels. The statistical noise is still large but it seems 
that the overall behavior does not depend too much 
on the kind of trial functions, although small de- 
tails could still be different. Note, however, that 
the liquid has little structure. More investigations 
of the influence of the trial function on the liquid 
structure is certainly needed, in particular, in the 
molecular dissociation region. 

4. Conclusions 

We have reported important progress in CEIMC, 
an efficient and accurate method to perform ab- 
initio simulations of condensed system with QMC 
energies. We have shown how the method performs 
in the case of hydrogen at high pressure, the sim- 
plest, but yet not understood, system. The new 
method allows us to cover a range of temperatures 
inaccessible to previous QMC methods for hydro- 
gen, a range where most of the interesting physics 
of hydrogen occurs, including the melting of the 
molecular and proton crystals, the molecular disso- 
ciation both in the liquid and in the crystal and the 
metallization of the system. 

A key ingredient in CEIMC is the trial function 
used to represent the electronic ground state. Even 
when a projection technique such as Reptation 
QMC is exploited to improve the bosonic part of 
the trial many-body wave function, its fermionic 
part, that is its nodal surface, is still playing a very 
crucial role in determining the electronic energies 
and therefore the overall thermal behavior of the 
system. In the present paper, we have reported a 
detailed investigation of these effects for hydrogen 
by comparing a number of different trial wave func- 
tions at two densities. We have shown as a fully 
analytical trial wave function, that is optimal in 
terms of computational efficiency in CEIMC, and 
which has been previously demonstrated to provide 
excellent accuracy for crystalline states, degrades as 
soon as some disorder is introduced in the protonic 
configurations. This result has been established 
by comparing with results for new trial functions 
obtained from a Slater determinant of IPP/LDA 



orbitals together with a two-body Jastrow correla- 
tion factor. A further backflow transformation of 
these orbitals has been introduced and character- 
ized. The new trial functions provide lower energies 
and more uniform overlap over a number of fixed 
representative configurations, which we use as an 
indication of the overall quality of the trial function. 
The most striking result on disordered configura- 
tions is that the LDABF energies at the VMC level 
are lower than the fully projected energies from 
the metallic trial function. This indicates that the 
improvement of performance comes mainly from 
the different nodal surfaces, while the bosonic part 
is responsible only for smaller improvements. The 
failure of the metallic wave function is most prob- 
ably due to the presence of some degeneracy of its 
orbital structure around the Fermi surface which is 
removed by solving the instantaneous band struc- 
ture. On the other hand, the use of complex wave 
functions and twist averaged boundary conditions 
in connection with the metallic trial function was 
expected to remove most of these degeneracies. A 
better understanding of this failure is desirable and 
deserves more investigation. 

The difference in energies for different trial func- 
tions, or more precisely the dispersions of the ener- 
gies from different wave functions, translates into a 
overall temperature factor at thermal equilibrium. 
The metallic trial function provides a dispersion 
which is roughly twice that of the corresponding dis- 
persion from the LDABF trial function. Therefore 
the metallic g PP {r) at temperature T should corre- 
spond to the LDABF g pp (r) at ~ T/2. This is indeed 
observed and the new g pp (r)'s from LDABF are in 
fair agreement with predictions of Car-Parrincllo 
MD|15j. This agreement remains somehow surpris- 
ing since, beyond the different methods of sampling 
protonic configurational space, the electronic de- 
scription in the two methods is quite different. We 
use LDA orbitals with a backflow transformation 
and a two body RPA Jastrow while in CPMD, only 
LDA orbitals are employed. Adding the backflow 
and the Jastrow, we obtain a fair gain of energy 
and, moreover, we can improve the bosonic part of 
the trial function by projecting in imaginary time. 
Further, we strongly reduce the finite size effects 
by averaging over the undetermined phase of the 
wave function, while CPMD calculations are per- 
formed at the r point only for closed shell systems 
(N p = 54 and 162). However the final agreement 
between the two methods indicates that the effects 
of these improvements on energy differences is only 
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minor. On the other hand, it is well known in simple 
liquids that g(r) is not very sensitive to changes of 
the interaction potential and this might explain the 
observed agreement. 

At lower densities, employing IPP orbitals and 
RPA Jastrow, we have recently found[14j a contin- 
uous molecular dissociation with density, at vari- 
ance with CPMD which has predicted a first order 
molecular dissociation transition[16j. The reliabil- 
ity of IPP trial function was only tested on crystal 
structures and should be further investigated for dis- 
ordered configurations along the lines shown here. 
This study is in progress. A recent BOMD study [IT] 
within DFT / GGA has reported a continuous molec- 
ular dissociation in agreement with our findings. 
This agreement suggests that improving the trial 
functions from IPP to LDABF might change the de- 
tails of the results but not the overall picture. This 
confirms that our present method can be most useful 
in condition where new interesting physics is hap- 
pening, such as near a liquid-liquid phase transitions 
or a metallization transition. 
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